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Abstract 

The hydrodynamic part of the velocity autocorrelation function of a granular fluid in the ho¬ 
mogeneous cooling state has been calculated by using mode-coupling theory for a finite system 
with periodic boundary conditions. The existence of the shearing instability, leading to a diver¬ 
gent behavior of the velocity flow fluctuations, is taken into account. A time region in which the 
velocity autocorrelation function exhibits a power law decay, when time is measured by the num¬ 
ber of collisions per particle, has been been identified. Also the explicit form of the exponential 
asymptotic long time decay has been obtained. The theoretical prediction for the power law decay 
is compared with molecular dynamics simulation results, and a good agreement is found, after 
taking into account finite size corrections. The effects of approaching the shearing instability are 
also explored. 
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I. INTRODUCTION 


Q 


Almost forty five years ago, Alder and Wainwright [l| reported, from a molecular dy¬ 
namics simulation study, the observation of an asymptotic power law decay of the velocity 
autocorrelation function (VACF) of a tagged particle in a fluid at equilibrium. At long times, 
the VACF decays as where d is the dimensionality of the system. They also proposed 

a simple hydrodynamic explanation, suggesting that the power law decay is due to the slow 
relaxation of the velocity vortex that is generated by the motion of the tagged particle in 
the fluid. Theoretical analysis of this effect have been carried out since then using kinetic 
theory 2, ^ and also more phenomenological mode coupling theories j4|. One of the main 
implications of the long time tails is that the time independent Navier-Stokes transport co¬ 
efficients, as dehned by the Green-Kubo relations, do not exist in two-dimensional systems, 
outside the limit of dilute gases, where the long time-power law tails disappear. 

Granular gases have attracted a lot of attention in recent years, not only because of the 
rich phenomenoly they exhibit, but also because they are considered as a proving ground 
for kinetic theory and non-equilibrium statistical mechanics. The methods used for elastic 
molecular fluids have been extended to the case of particles colliding inelastically j^, 1^ . In 


particular, hydrodynamic equations to Navier-Stokes order have been derivec 
Green-Kubo like expressions for the transport coefficients, both for dilute j?. 


granular gases j^, [ij. The low density expressions have been evaluated numerically by 


with formal 
8| and dense 


using the direct Monte Garlo simulation method [U, Q and a satisfactory agreement has 
been found with the results obtained by the Ghapmann-Enskog procedure applied to the 
inelastic Boltzmann equation in the hrst Sonine approximation [l3|. Note that in the low 
density limit considered in these works, algebraical decay of the correlation functions are not 
expected, by analogy with molecular gases. We are not aware of any numerical evaluation of 
the formal expressions of the transport coefficients for dense granular gases, aside from the 
self-diffusion coefficient in a three dimensional system Gonsequently, it is not known 
whether the Green-Kubo like expressions for the transport coefficients of granular gases do 
actually exist beyond the low density limit in two-dimensional systems. 

Due to dissipation in collisions, isolated granular gases (or with periodic boundary con¬ 
ditions) do not have an equilibrium Gibbs state, but rather there is a time-dependent ho¬ 
mogeneous cooling state (HGS). The Navier-Stokes transport coefficients of a granular gas 
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are expressed as averages over the distribution function corresponding to the HCS l7Hio|. 


& 


Consequently, it is a relevant question to know whether time correlation functions computed 
in this state exhibit slow decaying long time tails. 

Studying the existence of hydrodynamics in a two-dimensional molecular system can ap¬ 
pear as a rather formal and academic issue. Nevertheless, the situation is quite different 
when dealing with granular fluids. To reach and maintain them in a steady state, a perma¬ 
nent energy supply is necessary. The usual experimental procedures are either by means of 
an external held or by injecting energy through the boundaries. The price to be paid is that 
the system becomes highly inhomogeneous. Recently, a new possibility is being explored. 
The idea is to place a granular gas between two large parallel horizontal plates separated by 
a distance larger than one particle diameter, but smaller than two, in such a way that the 
system is actually a granular monolayer [^, 16|. To keep the system huidized the horizontal 
plates are vibrated. Then, the two-dimensional dynamics obtained by projecting the motion 
of the grains on an horizontal plane is observed. Experiments show that the behavior of the 
projected system resembles that of a (two-dimensional) fluid. Developing a self-consistent 
hydrodynamics for it seems an interesting and promising topic. Of course, this requires 
the investigation of the decay of the fluctuations and correlations in the two-dimensional 
system. Actually, an experimental setup very similar to the one described above has been al¬ 
ready used to measure velocity correlations on the hydrodynamic scale in a two-dimensional 
granular gas 


Long time tails in granular gases have already been studied. Kumaran 18|] considered 


sheared granular flows and reported the suppression of the power-law long time tail. Also, the 
tail of the velocity correlation function has been analyzed in the stationary state generated by 


submitting the grains to an external noise force or thermostat 


19| . Although the results are 


interesting, it is worth to remark that the possible relationship of this kind of forces with real 
experiments has not been established, and that the results depend strongly on the specihc 


form of the used force [19|. The case of an isolated granular gas has been investigated by 
Ahmad and Puri , by means of very extensive molecular dynamics simulations. Although 
their results for the velocity autocorrelation function in the two-dimensional case seem to 
suggest the existence of a power-law long time tail of exponent —1, when time is measured 


by the number of collisions per particle, they can not be considerec 
no comparison with theoretical predictions is carried out. In ref. 


as conclusive, and 


2l| . several current 
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autocorrelation functions have been investigated also in freely evolving granular fluids. Using 
a mode-coupling theory it is found, in particular, that the VACF has a long time decay 
of the same form as in molecular systems, when again time is measured by the number 
of collisions per particle. Moreover, extensive molecular dynamics simulation results are 
presented, and it is claimed that they are consistent with the theory, although again no 
quantitative comparison is carried out, other than the identification of a time region in 
which the behavior of the correlation function follows the power law predicted by the theory 
for the asymptotic long time limit. In the last two studies, very large systems are considered, 
actually inhnite in the theory developed in ref. j^, so that the HCS is highly unstable. In 
practice, this means that the state of the system departs from the HCS very soon, developing 


velocity vortices and later on strong density inhomogeneities |22|, 


23|. Consequently, the 


physical meaning of the VACF becomes rather uncertain after a very short period of time. 
In particular, there is no reason to expect it to be related in a simple way with the self¬ 
diffusion coefficient of the simulated system, if this coefficient exists. 

In this paper, the decay of the velocity autocorrelation function in a granular gas in 
the homogeneous cooling state (HCS) will be addressed. Both in the theory and in the 
simulations only systems in that state will be considered, although the effects of the shearing 
instability as it is approached will be taken into account. Therefore, in the systems analyzed 
here, the relationship between the VACF and the self-diffusion coefficient through a Green- 
Kubo formula is well established {q]. Although this implies the limitation to hnite systems, 
it is possible to extrapolate the results and identify scaled behaviors that are independent 
from the size. In particular, this happens with the existence of a time window for which 
the VACF has an algebraic decay that is correctly predicted by a mode-coupling theory, not 
only qualitatively but also quantitatively. 

The remaining of this paper is organized as follows. In Sec. m the dehnition of the 
self-diffusion coefficient in a granular gas in the HCS, and the steady-state representation 
of the latter, are summarized. It is important to stress that the existence of the steady 
representation is not an approximation, but an exact consequence of a change of variables. 
In particular, the self-diffusion coefficient is also determined by the VACF in the steady 
representation. In Sec. IIIIl the mode-coupling theory of Ernst, Hauge, and van Leeuwen 
jd] is rederived for a hnite granular gas with periodic boundary conditions in the HCS. The 
peculiarity of the hydrodynamic huctuations, playing a fundamental role in the theory, and 
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the effects of the shearing instability, are analyzed in Sec. IIVI The existence of a time scale 
over which the VACF has a power law decay in time is shown in Sec. |Vl where it is also 
discussed the exponential decay occurring in the asymptotic long time limit, due to the finite 
size of the system. In Sec. I VII the method used to take into account finite size effects when 
comparing molecular dynamics simulation results and theoretical predictions is described. 
The power law tails from the simulations are identified in Sec. IVIII Both, the power law and 
its amplitude are shown to be in quite good agreement with the mode-coupling predictions. 
Finally, Sec. IVIIII contains a short discussion of some relevant conceptual issues addressed 
in the paper, as well as some indications of possible extensions of the reported work. 


II. SELF-DIFFUSION COEFFICIENT IN THE STEADY-STATE REPRESENTA¬ 
TION OF THE HOMOGENEOUS COOLING STATE 

In Ref. j^, the self-diffusion coefficient D{t) of a granular gas of N inelastic hard spheres 
or disks of mass m in the HCS is identified from the diffusion equation, 

d 


dt 


ni(r, t) — D{t)W ni{r, t) = 0, 


( 1 ) 


where ni{r, t) is the number density of tagged particles at position r and time t. The formal 
expression of D{t) is given in terms of the VACF by the Green-Kubo formula 


1 


D{t) = ^ / dt' < vi{t) ■ vi{t');0 >HCS ■ 


( 2 ) 


Here, vi{t) is the velocity of a tagged particle at time t, and the angular brackets denote 
average with an initial distribution corresponding to the HCS of the system. Upon deriving 
the above expression, it is assumed that the system remains in the HCS in the time interval 
between 0 and t. Since all the particles of the system are mechanically equivalent, the 
VACF can be actually computed by using any of the N particles in the system. The 
HCS distribution function, Phcs, giving the probability density for finding the particles at 
positions {g*} with velocities {vi}, has the scaling form (2J, [25| 


pHCsiiqi} A^i} ,t) = [uoW] ^"^Phcs ’ 


where qij = Qi - Qj, 


Vo{t) = 


2T{t) 


m 


1/2 


( 3 ) 


( 4 ) 
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T{t) is the granular temperature (the Boltzmann constant is set equal to unity upon defining 
the granular temperature from the average kinetic energy), and p^cs ^ dimensionless 
isotropic distribution, which is invariant under space translations. In Eq. ([3]), all the time 
dependence due to collisional cooling occurs through the granular temperature, that obeys 
the Haff law 26] 

dvo(t) 1 , 

( 5 ) 

The cooling rate ({t) is proportional to vo{t). Then, it follows that the long time behavior 
of the thermal velocity vo{t) of a granular gas in the HCS is given by 


M't) ~ (C^) 


-1 


( 6 ) 


In this expression, a modified cooling rate coefficient. 


(7) 


2vo{t)' 

which does not depend on time, has been introduced. 

It is convenient to change to a new time scale in which the intrinsic time dependence 
of the HCS can be scaled out in some way. This allows to formally eliminate one of the 
two time explicit dependencies of the VACF in Eq. (|2]), making the theoretical analysis 
simpler and more direct 9|. Moreover, a direct molecular dynamics simulation of an freely 
evolving granular gas, as described by the dynamics in the actual phase space variables, has 
the limitation that the typical velocity of the particles becomes very small rather soon and, 
therefore, numerical inaccuracies become very large. This is a very serious difficulty when 
the interest is focused on the long time behavior of a property of the system, as it is the 
case here. To circumvent this problem, the dynamics of a system of inelastic hard spheres 
or disks in the HCS is exactly mapped onto the dynamics around a steady state 
This is done by defining a new time scale s by 

t 


28] 


UqS = In 


tn 


( 8 ) 


where to ^md Uq are two arbitrary constants. Consistently, the velocity cjj of a particle i in 
the new time scale is given by 


^i(t) = ^ = UJotVi(t). 


(9) 
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The particle dynamics consists now of an accelerated streaming between collisions, 


( 10 ) 

d 

= u:QU}i{s), ( 11 ) 

while the collision rnle in the new variables is the same as in the original ones, as a con¬ 
sequence of the linearity of the transformation given by Eq. ([9]) and the instantaneous 
character of collisions. 

The distribution function of the HCS in the transformed phase space reads 

PHCs{{qi] p*HCS ({gp}, • (12) 

In the long time limit, v^it) is given by Eq. ([6]) and, therefore, pucs becomes independent 
from the time s, and takes the stationary form 

?..({9.}.{‘^.})= (y) fc({9«}.{^}). (13) 

Let us define the temperature T of the HCS in the modified dynamics by 

d /jX/ \ nz 9 . ■. 

-nT{s)=<—up,s>, (14) 


where n is the total number of particles density of the granular gas, and the angular brackets 
denote average in the HCS in the new phase space. 




{^i])pHcs{{qi] ■ ( 15 ) 


It follows from Eq. (IT^ that, after some transient time interval, the system reaches a steady 
state with a temperature given by 




(16) 


This relationship provides a very efficient way to measure the cooling rate of a granular 
gas in the HCS by means of numerical particle simulations using the steady representation 


27|,l28|. 


The expression of the self-diffusion coefficient in the HCS, as given by Eq. ([2]), becomes 
simpler in the steady representation, specially when the initial time is chosen such that 
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the asymptotic steady state has already been reached. Then, the Green-Knbo form for the 
self-diffusion coefficient in the steady representation is trivially obtained, 

D{t) = i (W-\ ' f ds' < a;i(s') ■a;i(0) >,* , (17) 

u \ Tst / Jo 

where the ensemble average now is done with the stationary distribution reached with the 
modihed dynamics in the long time limit. In the above equation, the two time scales s and 
t are related by Eq. (|8]). 


III. HYDRODYNAMIC COMPONENT OF THE VACF 


Equation (HHi shows that the relevant VACF for the calculation of the self-diffusion 
coefficient in the steady representation is 

Cajuj{s) = ^ < ^l(s) ■ ^l(O) >st ■ (18) 

In this function, the trivial dependence on time occurring in Eq. ([2]) through the temperature 
of the HCS has been eliminated. In the following, attention will be focused in analyzing 
the long time behavior of Ci^^ois), that is important for the calculation of the self-diffusion 
coefficient as well as for the existence of the coefficient itself. 

The VACF is a spatially homogeneous function, but it can be expressed as the 

volume integral of an inhomogeneous quantity, 

Cu;Us) = ^ J dr j dr' < J{r,s) ■ J (r', 0) >su (19) 

where J is the microscopic current density of the tagged particle. 


J{r,s) = Wi5[r - qi{s)]. (20) 

Next, following the ideas of Ernst, Hauge, and van Leeuwen jd, 29|, the ensemble average 
in Eq. flT^ will be performed in two steps. First, a partial average for hxed values of 
the “relevant” magnitudes is taken. Afterwards, the average over the fluctuations of these 
magnitudes is carried out. The relevant magnitudes in the present context are those which 
are coupled to the velocity of the tagged particle and relaxing slowly. Here, it will be assumed 
that they are the same for the dynamics of a tagged particle in a granular gas in the HCS 
as for self-diffusion in an elastic molecular fluid at equilibrium, namely the local number 





density of the tagged particle and the local momentnm density. The reason is that they are 
conserved qnantities and, therefore, are expected to decay on a slow or hydrodynamic time 
scale. 

The microscopic density of the tagged particle at point r is given by 


Ni{r) = 6{r - qi). (21) 

Instead of the local momentnm density, the local velocity flow W{r) will be employed. It 
is dehned as 

1 ^ 

= (22) 

i=i 

In Fonrier space, the above two qnantities read 


N 


n 


ujje 


ikqj 


(23) 

(24) 


i=i 


A sqnare (d = 2) or cnbic (d = 3) system of side L will be considered, and periodic 
bonndary conditions employed. Now, a constrained distribntion in the modihed phase space 
T = {qi, (jJi, z = 1,..., iV}, for given macroscopic helds rii and u, of the relevant magnitndes, 
is dehned as 




Uk ^ - Nik) S {Uk - Wk) M^) 

P{{nik,Uk}) 


where 




a)S(u^-Wt)p„(T). 


(25) 


(26) 


The distribntion P{{nik, Uk}) can be nnderstood as the probability density of a hnctnation 
of both the nnmber of tagged particles density and the local velocity helds. Note that the 
constrained distribntion dehned by Eq. fl25l) is trivially normalized in the modihed phase 
space. 

The time dependence in the VACF can be made explicit by means of the psendo- 

Lionville operator C of the grannlar system in the scaled dynamics [9| 

1 


a^js) = 


d 


dr / dr' / dF J(r) • e 


Jir')Pst{r) 


(27) 
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Using the constrained distribution defined by Eq. fl2^ this can be rewritten as 


Cujoj{s) = - 


Y\_dnikduk P{{nik,Uk}) 


X dr dr' / dT J(r) ■ e 


{nik,'u.k}) 


( 28 ) 


Consider the average of the current density in the restricted steady HCS ensemble, 

/'dfj(r)^,f(f;{nifc,iifc}). (29) 


It seems sensible to assume that 


j8'>(r') = ni(r)u(r), 


(30) 


where ni{r) and u{r) are the inverse Fourier transformed of riik and u^., respectively. Then, 
using the Parceval relation, it is found 

J drj^'^^r) = ^^nikU^k, (31) 

k 

with V = L'^ the volume of the system. The crucial hypothesis of the theory will be intro¬ 
duced at this point. The right hand side of Eq. fl28|) is split into a fast, kinetic relaxation 
towards a local steady distribution followed by a much slower relaxation controlled by hy¬ 
drodynamics, 

CbQ;^kin(s) -f- (s). (32) 


The concept of local steady state is a direct extension of the widely used local equilibrium 
state of molecular systems. It is a reference state in which the system is considered to be 
locally in the steady HCS, but with the hydrodynamic fields density of tagged particle and 


flow velocity, being functions of space and time 


30| . In the following, attention will be 


restricted to the second term on the right hand side of the above equation, that is expected 
to dominate for s S> Srei, where s^ei is some characteristic microscopic relaxation time, i.e. 


— ^ujcj,hyd(^) j 


(33) 


for s S> Srei- Moreover, in the spirit of the above time scale separation, it is assumed that 
in the hydrodynamic regime. 




~ j^^^(r')e {nik,Uk}). 


(c)l 


— sC~:{c) , 


(34) 
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Substitution of this expression into Eq. (j28|) . and use of Eq. (j^ . yields 
C^tA;a;,hyd ('5) — ?/ n dnikduk I -P({nifc, Uk}) 

k± \ k / 

xnikiU_ki- j dr j (irj(r)e"*'^p^t^(f ; {nifc, -Ufc})- (35) 

The next approximation is based once again on the same physical picture. Since the density 
of the tagged particle and the local velocity flow evolve very slowly for small values of the 
wavenumber vector, the distribution function 'pf^ adjusts itself continuously, remaining with 
the same functional form on the long time scale, 

; {nifc(s), Ufc(s)}). (36) 


Moreover, nik{s) and Uk{s) are determined by the linearized hydrodynamic equations (to 
Navier-Stokes order) with the initial conditions riik and Uk- This seems legitimate for small 
enough wave-vectors k. Thus it must be verihed a posteriori whether the long time behavior 
of is actually governed by small wave numbers. When Eq. fISBD is substituted into Eq. 
fl35|) one gets 

C'tjaj,hyd(’S) — n dfiikdukj P(^{nik, Uk}')TiikiU—ki j drj^^\r,s), ( 37 ) 

ki \ k / 

where Eq. fl2^ has been employed. If now Eq. fl3l|) is taken into account, the above expres¬ 
sion can be rewritten as 

C^ui,hyd{s) - J iY\_dnikdUk\pi{nik,Uk})niki-u_kinik2is)u-k2{s)- (38) 

fci fc2 \ k / 

In linear hydrodynamics around the HCS, the fluctuations of the density of the tagged 
particle and of the velocity flow are not coupled, i.e. they are statistically independent. 
Moreover, different wave-vectors k are also uncoupled, so that 


P{{nik,Uk}) = YlPniinik)Pu^iuk), (39) 

k 

where the marginal probability distributions P^ and verify the normalization conditions 

/ dnikPniinik) = / dukPu^,{uk) = 1. (40) 


In addition, the isotropy of fluctuations in the HCS implies that the only non-vanishing 
contributions in Eq. fl55]) are those with ki = —k 2 . In this way, it is obtained 


c'^^,hyd(s) 


dUkPni{nik)Pu^{uk)nikni_k{s)uk ■ u_k{s). 


(41) 
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IV. HYDRODYNAMIC FLUCTUATIONS 


As indicated above, the linearized hydrodynamic equations for a granular fluid will be 
used to evaluate the time evolution of the fluctuations of the hydrodynamic helds. The 
diffusion equation ([T]) in the s-time scale has the form 

d 


+ Dk^ nifc(s) = 0, 


with 


ds 


D= dsCi_ja}is). 

Jo 


(42) 


(43) 


Here it has been assumed that the VACF decays fast enough so the time integral on the right 
hand side of the above equation exists, leading to a time-independent transport coefficient 
D. Integration of Eq. (H2l) gives 


ni_fc(s) = ni_fc(0)e 


-k'^Ds 


(44) 


The flow velocity fluctuations Up, can be decomposed into their longitudinal and transversal 
components, Uk\\ and Uk±, dehned by 

Uk ■ k 


'^kW 


'^kl. — '^k 


Uk ■ k 




k. 


respectively. In the new time scale s, u± obeys the closed equation 


d 


13l 


■^-uo + vk'^] Wfc±(s) = 0, 


(45) 

(46) 

(47) 


where rj is dehned from the shear viscosity r] by 


V = 


T] 


mn 


Tst ^ 

nt)} 


1/2 


The solution of Eq. flT7)l is 


'^fc±('S) 


{uJo—rik?)s 


(48) 


(49) 


In this expression, it is manifest that perturbations of the modihed transversal velocity grow 
in time for small wave vectors, i.e. in large enough systems. This is the origin of the shearing 


instability of the HCS 


22 


23| . in which strong density inhomogeneities are developed. 
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They are generated by nonlinear coupling contributions of the transversal velocity 311 • To 
avoid misunderstandings, it is important to realize that the linear stability criterion for Uk,± 
following from Eq. 0491) does not depend on the arbitrary value of ujq. The critical value of 
the wavenumber vector /cx, such that transversal modes with /c < fcx are unstable, is given 
by the solution of the equation 

ujo — rjk^ = 0, (50) 

that, using Eq. (IT6|) . is seen to be equivalent to 


c- 




1/2 


2mT{t)) 


rjk^ = 0 . 


(51) 


It is now clear that k± does not depend on uq or T(t), since rj is proportional to T(t)^/^. 

Now, the contributions to C^a;,hyd(’S) from the fluctuations of the longitudinal component 
of the flow held Uk\\ should be considered. Determining its time evolution using the linear 
hydrodynamic equations for a granular huid is a rather involved problem, since Uk\\ can not 
be associated to a unique hydrodynamic mode, and its functional form changes depending 
on the range of values of the wavenumber considered ISj. Here the focus will be on small 
systems, in which the allowed values of A;, i.e. compatible with the periodic boundary 
conditions, are such that the system exhibits a sound-like hydrodynamic mode, at least in 
dilute systems 0. Then, .t wn, be assumed that the eonttibntion ftotn the longitudinal 
component of the how veloci^ to C^joj^s) for large s is sub-dominant, as it is the case in 


systems of elastic particles jd, 291. Alternatively, neglecting the contributions from the 


longitudinal component of u, can be considered as an incompressible huid approximation. 
Substitution of Eqs. (Hij) and into Eq. flTT|) and use of 


dnik P{nik)\nik\‘^ = 1 , 


leads to 


where 




a;a;,hyd 


5 ~ 




Vd 




(52) 


(53) 


A{k) = ^ j dukP{uk) \uk\^. (54) 

In Eq. fl53l) the range of values of k is restricted hj km < k < kui where km is set by the 
system size, km = ‘2tt/L, and kM is of the order of 27r times the inverse of the mean free 
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path. This upper bound condition is needed to guarantee that the local HCS varies slowly 
in space. 

In elastic, molecular fluids in thermal equilibrium, the calculation of the second moment 
of the velocity held fluctuations is trivial. However, in a granular gas in the HCS, velocity 
correlations are present, rendering the computation of A{k) more difficult. Using inelastic 
fluctuating hydrodynamics, van Noije et al. |32| computed the quantity 

1 


SAkA) = 


< Uk±{t) ■ Uk±{t)-, 0 >HCS ■ 


(55) 


V{d-1) 

Here u is the local flow velocity fluctuation in the original time scale t. The result, in the 
notation used in this paper is 

„2uJos(l-k^ fk^) 


S±{k,t) = 


T{t) 


mn 


1 + 


1 - kyk] 


(56) 


An equivalent result was obtained by using a single relaxation model kinetic theory 33|. 


Moreover, in recent years, a description of fluctuating hydro dynamics in dilute granular 
gases more rigorous and complete than the one used in ref. |32j has been developed [34 1- 
Nevertheless, for the values of the coefficient of normal restitution to be considered here, 
that correspond to the quasi-elastic limit, both theories lead to results quantitatively indis¬ 
tinguishable for S±{k,t)- For systems in which the HCS is stable, it is k^ > k±, and the 
long time limit of the above expression is 

T{t) e 


S±{k,t) = 


k'i' 


(57) 


mn k'^ 

It is worth to point out that the result given by this equation has the scaling property 
implied by the assumed form of the distribution function of the HCS, Eq. ([3]), while Eq. (]56|1 
does not. Then, taking into account the relationship between the velocity fields in the time 
scales t and s, it is identified 

{d - i)fst e 


A{k) = 


mn 


k? 


q' 


( 68 ) 


Introduction of this into Eq. fl5^ provides the explicit expression for the hydrodynamic part 
of the steady VACF of a granular gas in the HCS, 


Cajaj,hyd ('5) — 


mnV d 






(59) 
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The prime in the sum over k indicates the restriction imposed by the two cutoffs mentioned 
before. This expression leads to an exponentially increasing VACF on the time scale s for 
values oi k < kc = oJo/{ri + D). Nevertheless, this is not physically relevant since the HCS 
is unstable in systems large enough as to allow wave-vectors with these values. This is 
because kc < k±, where k±, dehned by Eq. flSUD . is the lower bound of the wave-vector for 
the stability of the HCS with regards to the shearing instability. Moreover, in Eq. (I59|) 
Ca;aj,hyd(s) also exhibits a divergent amplitude when km approaches k±. This divergence is 
just the shearing instability, and clearly implies a divergent behavior of the self-diffusion 
coefficient as given by the Green-Kubo expression, Eq. ([2]), as the instability is approached. 


V. TWO LIMITING TIME REGIMES OF THE VACF 


In this section, two particular limits of Eq. fl59|) will be investigated. Suppose hrst that 
L -C Lc, where is the critical size of the system for the shearing instability, i.e. (see Eq. 

L, = 2w ’’ , (60) 

/ ~ , \V2 

where Vst = ( 2Tst/mj . Then it is. 




27r 


> 1, 


for all the allowed values of k. It follows that Eq. (1591) can be approximated by 

{d - v-W 




a;a;,hyd 


5 ~ 




(61) 


(62) 


mnVd 

A useful representation of this expression can be obtained by means of the d-dimensional 
Poisson-sum formula 


35| 


^ I 


—2'KiLl-r 


S(r), 


(63) 


where n and Z are d-dimensional vectors whose components are integers, the summations 
extend from —cx) to -|-cx) for each of the components of n and /, and the r integration 
extends over the inhnite d-dimensional space. Indeed, use of Eq. fl63|l into Eq. fl62|) yields 

d/2 

e 4(f7+o)s _ (^ 54 ) 
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Consider times s such that 

s < So = —7--. (65) 

4 f?7 + Dj TT^ 

Taking into account that its has been supposed that L Lc and Eq. (I50|) . it is easy to 
verify that the above condition also implies that 

ojqS C 1 . ( 66 ) 


As a consequence, Eq. 


reduces to 


("S) — 


{d - 1)T,, 
mnd 


4:71 (7] + D) s 


d/2 


(67) 


This prediction is expected to be valid for short times on the hydrodynamic scale, in the 
sense of being s <C Sq, but, on the other hand, the equation only holds after all the fast 
non-hydrodynamic modes have decayed. It is important to compare this result with the 
long time tails the of the VACF in an e^ilibrium molecular system. Formally, the derived 
expressions in both cases look the same |4, l29|, but there are relevant conceptual differences. 
The time scale s used in Eq. (167)) is related with the original time scale t by Eq. ([8]), so 
that the algebraic decay transforms into an even slower logarithmic decay on the time 
scale t. Actually, the time scale s is proportional to the cumulated number of collisions per 
particle for a system of inelastic hard spheres or disks in the HCS {q]. In a molecular system 
in equilibrium, this number is just proportional to t. Another, more signihcant, difference 
is that Eq. (167)) has been obtained here as valid in an intermediate time regime, while 
the corresponding expression in a molecular system has been usually derived as valid for 
asymptotically long times or, more precisely, for all the hydrodynamic decay of the VACF. 
The reason for this strong difference is that in molecular systems, the thermodynamic limit 
was considered, while in the present case such a limit can not be taken, at least in the usual 
way, for the HCS due to the shearing instability. A particularly enlightening discussion of 
the hnite size effects in the evaluation of the VACF in a molecular system can be found in 
|. A hnal comment on Eq. fl67)) seems appropriate. Although the time scale in which 


ref. 


the behavior predicted by this equation holds depends on the size L of the system, the shape 
of the VACF in that time region is independent of L, as long as it is well inside the range 
in which the HCS is stable with regards to the shearing instability. 
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The second relevant limit of Eq. fl59|) to be considered, is its asymptotic behavior for very 
large times. Then, it will be supposed now that s S> so, with sq given by Eq. fl6^ . In this 
limit, it is evident that contributions with the smallest possible value of k, km, will dominate 
in Eq. flS^ . Taking into account that, in a square or cubic geometry with periodic boundary 
conditions there are d modes with k = km, one gets 


C*a;t<;,hyd ('5) — 


mnV km — k^^ 


( 68 ) 


Note that the stability condition for the HCS, rjk'^ > ojq, implies that always 

decays exponentially in the asymptotic long time limit in the scale s. On the other hand, 
the amplitude of the decay diverges as the instability is approached. 


VI. MOLECULAR DYNAMICS SIMULATION RESULTS FOR THE VACF 


In order to check the accuracy of the theory developed in the previous sections, event 
driven molecular dynamics (MD) simulations of a system of inelastic hard disks have been 
performed. The coefficient of restitution has been always taken as a = 0.99, while the 
number density has been varied in the interval 0.231 < ncx^ < 0.385. The maximum 
and minimum density values in the interval correspond to densities studied by Alder and 
Wainwright in their seminal paper [l|, so it is possible a direct comparison with their results 
for elastic systems. Simulations with different numbers of particles in the range 500 < N < 
1000 have been carried out for each density, modifying accordingly the size L of the system. 
By estimating the critical wave number for the shearing instability of the HCS, k±, using 
the transport coefficients and the cooling rate obtained by means of the revised Enskog 
theory for inelastic hard disks 37|, it is found that in all the simulations considered, the 
minimum wave number allowed, km, is significantly larger than k±. Moreover, it was always 
checked in the simulations that the system remained homogeneous and with no appreciable 
velocity vortices. This latter condition is specially relevant in the present context, since the 
fluctuations of the transversal velocity play a dominant role in the long time, hydrodynamic 
behavior of the VACF. Attention has also been paid to consider times such that the particles 
do not interact with their images in other cells. To avoid it, Erpenbeck and Wood 36 1 


suggested to consider times smaller than the time a sound wave takes to travel a distance 
equal to the linear size of the system L. This criterium is adopted in most studies of the 
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VACF, see for instance |38j, and has also been nsed here. 


In the simnlations, the steady representation of the HCS discnssed in Sec. [TTl has been 
nsed. Althongh it is not relevant for the resnlts presented below, let ns mention that the 
valne of the parameter uq was chosen snch that, if the Enskog theory were exact, the valne 
of the steady temperatnre, as predicted by Eq. (fTBD . wonld be the same as the initial one 
Tst = T{0). In all the simnlations, the applied procednre was as follows. The system was 
prepared in an initial spatially homogeneons state, with a Ganssian velocity distribntion. 
This initial state was allowed to evolve with the modihed dynamics and periodic bonndary 
conditions, nntil a steady state was reached. Then, the VACF was measnred. The resnlts 
were averaged over several trajectories, typically a few hnndreds. 

As already pointed ont, the prediction abont the existence of a hydrodynamic period over 
which the VACF exhibits a power law decay in the s time scale, and the specihc form of 
this decay, do not involve the size of the system, aside from repairing it to be smaller than 
the critical size for the shearing instability, as it can be seen in Eq. fl67ll . Nevertheless, when 
compnting the VACF from MD simnlations, there are two main (related) reasons for which 
the hnite size of the system mnst be taken into acconnt to compare nnmerical resnlts with 
the theoretical predictions. The hrst reason is dne to conservation of the total momentnm 
of the system 39|, that implies that a given particle is actnally moving in a flnid with a 
non-zero average velocity. The second canse is the difference between the ensemble in which 
the theory is developed and the one nsed in the simnlations 36|. Both effects have been 
analyzed and qnantihed for molecnlar systems in eqnilibrinm, bnt the argnments nsed in this 
case are not easy to extend to the present intrinsic non-eqnilibrinm system with the particles 
snbmitted to an effective acceleration. Then, an henristic view will be adopted. It will be 
assnmed, as it is the case in elastic systems, that both effects scale with N~^ to leading order, 
and the proportionality constant will be determined from the simulations themselves. Of 
course, the simulations also provide a test, a posteriori, of whether the leading dependence 
on N is indeed the assumed one. 

In Fig. [H the results obtained for the decay of the VACF in a system with density 
= 0.385 and different values of the number of particles N are shown. The quantity 
actually plotted is the dimensionless VACF, fui{s), dehned as 

mC^Js) 


UM = 


(69) 
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FIG. 1; (Color online) Dimensionless normalized VACF fu}{s) for a system of inelastic hard disks 
in the HCS. Time s is measured in units of (T{0)/ma‘^^ , where r(0) is the initial granular 

temperature of the system. The coefficient of normal restitution is a = 0.99 and the density 
is na'^ = 0.385. Different numbers of particles (and sizes of the system) have been used in the 
simulations, as indicated in the inset. 

Although small, the dependence of the results on the number of particles is clearly identified 
in the figure and, at a given time s, the VACF is smaller the smaller the system. To 
investigate the dependence on N, three different times, representative of the hydrodynamic 
relaxation of the VACF, were considered. At each of these times, the values of the VACF 
obtained in systems with different sizes but the same density, were analyzed. The results 
are given in Fig. |2l The solid lines are linear fits at each time, the slope of the fitting being 
very similar for the three times considered, namely 2.93 ± 0.08. This slope is about three 
times the one predicted and observed in elastic systems at equilibrium, showing that the 
dependence on the number of particles is larger for inelastic particles in the HCS than for 
molecular systems at equilibrium 

Denoting the VACF measured in the MD simulations with N particles by foj,N{s), h is 
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FIG. 2: (Color online) Values of the dimensionless normalized VACF as a function of the inverse 
of the number of particles N of the system. The coefficient of normal restitution is a = 0.99 and the 
density tkt^ = 0.385. The symbols are simulation results at three different times, namely s = 2.6 
(black), s = 5.0 (red), and s = 7.2 (green), from top to bottom. Time s is measured in units of 
ma^ 1 , where T(0) is the initial granular temperature of the system. The straight lines 

are linear fits of the simulation data. The values of the slopes of the three lines are very close. 


concluded that the extrapolated value, for an infinite system, is given by 

Us)^Un{s) + 2.93N-\ (70) 

The modified curves, obtained by adding this correction to the measured VACF reported in 
Fig. [T] are plotted in Fig. [31 It is observed that the collapse of the several curves is very 
good, especially for s < 20 iT{0)/ma‘^] . Similar results have been obtained for other 

densities in the interval considered. Moreover it is found that the prefactor of N~^ in the 
finite size correction term increases as the density increases. On the other hand, in the 
elastic case, it is always equal to unity, independently of the density. 

Before proceeding any further, the above limiting process must be put in a proper con- 
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FIG. 3: (Color online) The same as in Fig. [H but now the finite size correction given in Eq. (j70[) 
has been added to each curve. 

text. It is not claimed that the derived results hold in the limit of an inhnite system. This is 
not true, since an inhnite system of inelastic hard spheres or disks can not exist in the HCS, 
because it is very unstable. What has been done is to take into account the conceptual dif¬ 
ferences between theory and MD simulations, due to the use of different physical conditions 
as indicated above, realize that they scale with the inverse of the number of particles, and 
eliminate those differences to carry out a fair comparison. 


VII. THE POWER LAW TAIL OF THE VACF 

Once it is known how to translate the MD numerical results for the VACF into results 
corresponding to the conditions under which the theory is developed, the existence of a time 
region in which the VACF presents a power law time tail, as predicted by Eq. fE7)) . will be 
investigated. As mentioned in the previous section, the observation time in the simulations 
is limited by the size of the system. Since, once corrected, the resulting VACF does not 
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depend on the number of particles used, only results with the largest number of particles 
compatible with the stability of the system, namely N = 1000, will be presented from now 
on. In any case, only results for times s shorter than the time it takes a sound wave to cross 
the system will be shown. 

In Fig. m the scaled VACF of a system of inelastic hard disks is shown, both as a 
function of time s/se and as a function of {s/se)~^, for two different densities, ncr^ = 0.231 
and ncr^ = 0.385, respectively. Here 


m \ 1 

TifstJ 2anx 

is the mean free time between collisions in the steady representation of the HCS and com¬ 
puted using the Esnskog theory, so that y is the equilibrium pair distribution for hard disks 
at contact. This is a convenient dimensionless time scale to compare results corresponding to 
different densities, since it is proportional to the cumulative number of collisions per particle 
which is the relevant time scale for the relaxation of the system. As expected, the velocity 
correlations are more persistent, i.e. they decay slower, the denser the system. Moreover, a 
clear region exhibiting a linear dependence on is identihed, in qualitative agreement with 
the theoretical prediction. To carry out a quantitative comparison between the simulation 
results and Eq. fl67|l . the latter (with d = 2) is substituted into Eq. fl69l) . and the result is 
rewritten in the form ^ 

fuj,hyd — CX.D ^ ) (”^ 2 ) 

where the amplitude od of the tail is given by 




— -wy-:3Y • ('^3) 

Sttuse (h + ) 

Equation d72|) is of the form considered in the original paper by Alder and Wainwright [l|, 
and also by other authors j^. 

The slope of the hts to straight lines of the linear in {s/se)~^ regions observed in the 
decay of the dimensionless VACF foj{s) obtained in the MD simulations, provides the values 
of the tail amplitude. In Fig. [5l a/) is plotted as a function of the number density na^ 
of the system. The symbols are the results from the MD simulations. In addition to the 
results for inelastic disks with a coefficient of normal restitution a = 0.99 (red squares), 
MD values for a system of elastic hard disks (black circles) have also been included. The 
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FIG. 4: (Color online) VACF for an inelastic system of hard disks with a coefficient of normal 
restitution a = 0.99 and two different values of the density, as indicated in the inset. The results 
have been obtained in a system with 1000 particles, and finite size effects have been corrected as 
discussed in the main text. 
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latter agree with those reported by Alder and Wainwright many years ago [l|, l36j . The solid 


fl 


line is the theoretical prediction given by Eq. fl73|l . nsing the expressions for the inelastic 
transport coefficients derived in the Enskog approximation j^, and the Henderson valne 


4l| for the pair distribntion fnnction at contact. The consistency of nsing the bare transport 


coefficients to compnte the algebraical tails will be discnssed in the last Section of the paper. 
Eor the sake of completeness, these expressions are reprodnced in the Appendix. Actnally, 
the theoretical prediction for the amplitnde of the tail for an inelastic system with a = 0.99 is 
indistingnishable from the theoretical prediction for a system of elastic hard disks. However, 
the different physical natnre of the states considered in both cases mnst be kept in mind. 
In particnlar, it is worth to insist on that a decay of the scaled VACF corresponds to a 


decay {t In t) ^ of the VACF in the original time (and velocities) scale. As already known 




the agreement between theory and simnlations is qnite good in the elastic case. For inelastic 
systems, the comparison can still be considered as satisfactory, althongh it is clear that the 
effects of the inelasticity are mnch larger than predicted by the theory. The inelastic tail 
amplitnde is np to 15% larger than the elastic valne at the highest density considered. 

VIII. CONCLUDING REMARKS 


Using mode-conpling theory, an approximate expression for the hydrodynamic, long time, 
part of the relaxation of the VACF of a hnite grannlar gas in the HCS has been derived. 
The analysis parallels in several aspects the one carried ont for elastic systems, bnt there 
are many relevant qnantitative and conceptnal differences. First, the correlation fnnction 
considered here corresponds to the modihed dynamics dehned by Eqs. flTOD and (ITT]) , in which 
the particles are snbmitted to an acceleration, as a conseqnence of a change in the original 


time scale. The new time scale measnres the average nnmber of collisions per partic’ 
plays a crncial role for the stndy of response fnnctions in non-eqnilibrinm systems j?. 


e, and 

im 


The self-diffnsion transport coefficient of the grannlar gas can be expressed in terms of a 
one-time VACF compnted in a steady state only in the modihed dynamics. Then, it is 
the time behavior of the correlation fnnction in the modihed dynamics the one relevant for 
compnting the transport coefficient by means of the Green-Knbo expression. 

Another important diherence between molecnlar systems at eqnilibrinm and grannlar 
gases in the HCS lies in the meaning of taking the thermodynamic limit in the MD simn- 
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FIG. 5: (Color online) Tail amplitude as a function of the number density of the system. The 
symbols are simulation results: the (black) circles are for a system of elastic hard disks and the 
(red) squares for a system of inelastic hard disks with a coefficient of normal restitution a = 0.99. 
The solid line is the theoretical prediction using the mode coupling theory developed in the main 
text and the transport coefficients obtained from the revised Enskog theory. On the scale of the 
figure, the predictions for the elastic and the inelastic systems are undistinguishable. 

lations. At a theoretical level, the HCS can be considered and analyzed at arbitrary size, 
density, and valne of the coefficient of normal restitntion, and this legitimates its nse as a 
reference homogeneons state. On the other hand, in the simnlations a very serious limitation 
shows up due to the shearing instability. At hxed density, as the coefficient of normal resti¬ 
tution decreases, the size of the system for which the HCS remains stable decreases rather 
fast. At moderate densities, only granular gases of inelastic hard spheres or disks with a very 
small number of particles can be simulated in the HCS. This limitation is specially severe 
when the interest is on the long time behavior of a given property, as it is the case here. 
Actually, it is quite difficult to reach the time window for which the velocity-autocorrelation 
function has a power law decay for values of the restitution coefficient not very close to one. 
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To illustrate the effect of approaching the instability, in Fig. |6] the VACF of a system 
in the HCS with a = 0.98 is shown at three different times as a function of the inverse of 
the number of particles N. The density of the system is in all cases ncr^ = 0.3. For small 
enough systems {N < 700), a linear in N~^ fit is accurate, similarly to what happens for 
a = 0, 99 and = 0.385 in Fig. [2J The measured slope now is 5.17 ± 0.69, which is two 
times the value for a = 0.99 and the same density. On the other hand, for N > 700, a 
departure from the linear dependence is clearly observed. We expect this behavior to be due 
to the proximity of the instability and the presence of large fluctuations of the transversal 
velocity field, as predicted by the factor of in the addends in Eq. fl5^ . All 

the results reported in the figure correspond to states in which all the hydrodynamic fields 
were observed to stay homogeneous. Therefore, it seems that the hydrodynamic part of the 
VACF has a divergent behavior as the shearing instability is approached. The same must 
happen with the (apparent) self-diffusion coefficient. 

In two-dimensional molecular systems, it was realized that the analysis similar to the 
one reported in this paper, is internally inconsistent. The reason is that it is assumed 
in the derivation that a finite self-diffusion constant D exists. But the asymptotic t~^ 
long-time tail of the VACF implies that the Green-Kubo expression for the self-diffusion 
coefficient diverges as In t. It is worth to mention that more refined mode-coupling theories 
using a time-dependent expression for the self-diffusion coefficient have led to the prediction 
that the t~^ decay corresponds to intermediate times, while a slightly faster decay, namely 


as 


-1 


tVhit] is expected at 


numerical simulations 
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ater times. Nevertheless, this has never been confirmed by 
. The view adopted here, and consistent with the results 


reported, is that the influence of the time tail on the observed value of the self-diffusion 
coefficient remains negligible over a time scale going well into the time scale in which the 
tails can be observed. On this scale, D, given by the Green-Kubo expression, appears 
as constant, the contribution from the tails remaining very small. What happens at the 
far end of the region in which the self-diffusion coefficient seems to be constant, remains an 
open question, except by the mentioned theoretical predictions. 

A completely different question is the influence of the clustering instability on the observed 
self-diffusion coefficient in a finite system as the shearing instability is approached. This is 
a very interesting issue that deserves further investigation. 
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FIG. 6: (Color online) Values of the dimensionless normalized VACF /(^ as a function of the inverse 
of the number of particles N of the system. The coefficient of normal restitution is a = 0.98 
and the density mr^ = 0.3. The symbols are simulation results at three different times, namely 
s = 2.4 (black), s = 5 (red), and 10 (green), from top to bottom. Time s is measured in units of 
, where T{0) is the initial granular temperature of the system. The straight lines 
are linear fits of the simulation data in the region where this behavior is observed. The values of 
the slopes of the three lines are very close. 

IX. ACKNOWLEDGMENTS 

This research was supported by the Ministerio de Economia y Competitividad (Spain) 
through Grant No. FIS2014-53808-P (partially financed by FEDER funds). 

Appendix A: The cooling rate and the transport coefficients in the Enskog theory 


In this Appendix, the expressions for the cooling rate and the transport coefficients of 
a system of inelastic hard spheres or disks of mass m, diameter a and constant coefficient 


of restitution a, obtained by using the Enskog approximation 


37l | are given for the sake of 


27 






completeness. The cooling rate ( is 


c = c 


nT 

Vo 


where r/o is the elastic value of the shear viscosity in the dilute limit, 


Vo = 


{d + 2) T{d/2) 1/2 -(d-i) 

87r{d-i)/2 


The reduced cooling rate (* is given by 

with 02 the hrst coefficient of the Sonine expansion of the HCS distribution 

16(l-a)(l-2a2) 


24 


4^, 


0-2 — 


9 + 24d + (8d - 41)a + 30(1 - a)a^ 


The shear viscosity rj is 


T] = T] 


^d—1 


1 ++“> 


d 


where is the kinetic contribution to the viscosity, 

. 1 


V = Vo- 


2^d—2 

1-:—r(l + a)(l — 3a)0x 


K - K* L ^ + 2 

while 7 is the bulk viscosity, that vanishes in the dilute limit. 


7 = ho 


2^2d-\-l 
{d + 2)7r 


0\(1 + a) 1 
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In the above expressions, u* is 






and 0 is the volume fraction. 


TT 


d/2 


2'^r (1 + f) 




Finally, the bare self-diffusion coefficient is given by 451 

4 T{d/2)d 


D = 


rp \ 1/2 


(Al) 


(A2) 


(A3) 


(A4) 


(A5) 


(A6) 


(A7) 


(AS) 


(A9) 


( 10 ) 
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In the particular case of hard disks, d = 2, the pair distribution function at contact x is 


m 


X = 




1 - 
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